# ------------------------------------------------------------------------------#
# Description: Create Table A.3: Peer effects heterogeneity by visibility proxies
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(data.table)
library(fixest)
library(modelsummary)

options(scipen = 999)

# Load and merge base data -----------------------------------------------------
load("data/adoptions/adoptions.RData")
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale    := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale     := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale     := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale   := peer.adopt.both.100 / 100]

dat <- dat[year > 2010]
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Merge GSV visibility and sign data -------------------------------------------
load("data/gsv/rw_gsv.RData")

load("data/gsv/peer_adopt_buckets_any_vis.RData")
dat <- merge(dat, peer.adopt.bucket.any.vis, by = c("pin", "year"))
rm(peer.adopt.bucket.any.vis)

load("data/gsv/peer_adopt_buckets_any_rg_vis.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg.vis, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg.vis)

load("data/gsv/peer_adopt_buckets_cs_vis.RData")
dat <- merge(dat, peer.adopt.bucket.cs.vis, by = c("pin", "year"))
rm(peer.adopt.bucket.cs.vis)

load("data/gsv/peer_adopt_buckets_any_sign.RData")
dat <- merge(dat, peer.adopt.bucket.any.sign, by = c("pin", "year"))
rm(peer.adopt.bucket.any.sign)

load("data/gsv/peer_adopt_buckets_any_rg_sign.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg.sign, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg.sign)

load("data/gsv/peer_adopt_buckets_cs_sign.RData")
dat <- merge(dat, peer.adopt.bucket.cs.sign, by = c("pin", "year"))
rm(peer.adopt.bucket.cs.sign)

# Merge corner data ------------------------------------------------------------
load("data/corner/peer_adopt_buckets_any_corner.RData")
dat <- merge(dat, peer.adopt.bucket.any.corner, by = c("pin", "year"))
rm(peer.adopt.bucket.any.corner)

load("data/corner/peer_elig_bucket_cs_corner.RData")
dat <- merge(dat, peer.buckets.cs.corner, by = c("pin", "year"))
rm(peer.buckets.cs.corner)

load("data/corner/peer_elig_bucket_rg_corner.RData")
dat <- merge(dat, peer.bucket.rg.corner, by = c("pin", "year"))
rm(peer.bucket.rg.corner)

load("data/corner/total_corner_peers.RData")
dat <- merge(dat, total.corner.peers, by = "pin")
rm(total.corner.peers)

dat[, corner.peers.avg.100 := (corner.peers.cs.100 + corner.peers.rg.100) / 2]

# Merge bus data ----------------------------------------------------------------
load("data/bus/peer_adopt_buckets_any_bus.RData")
dat <- merge(dat, peer.adopt.bucket.any.bus, by = c("pin", "year"))
rm(peer.adopt.bucket.any.bus)

load("data/bus/peer_elig_buckets_cs_bus.RData")
dat <- merge(dat, peer.buckets.cs.bus, by = c("pin", "year"))
rm(peer.buckets.cs.bus)

load("data/bus/peer_elig_buckets_rg_bus.RData")
dat <- merge(dat, peer.bucket.rg.bus, by = c("pin", "year"))
rm(peer.bucket.rg.bus)

load("data/bus/total_bus_peers.RData")
dat <- merge(dat, total.bus.peers, by = "pin")
rm(total.bus.peers)

dat[, bus.peers.avg.100 := (bus.peers.cs.100 + bus.peers.rg.100) / 2]

# Merge census data -------------------------------------------------------------
load("data/census/census_bg.RData")
census.2014[, blkgrp := as.numeric(paste0(tract, bg))]
census.2014[, active.trans := pub.trans + bike + walk]
census.2014[, degree := college + adv.degree]

dat <- merge(dat, census.2014[, .(blkgrp, med.inc, non.white, non.white.asian,
                                  active.trans, degree, own.occ, gov.assist)],
             by = "blkgrp")
rm(census.2014)

dat[, hi.inc := fifelse(med.inc > median(dat$med.inc), 1, 0)]
dat[, hi.act := fifelse(active.trans > median(dat$active.trans), 1, 0)]
dat[, hi.own := fifelse(own.occ > median(dat$own.occ), 1, 0)]

# Estimate models ---------------------------------------------------------------
ivbase <- feols(adopt.any * 100 ~ peers.100 | blkgrp + year |
                  peer.adopt.any.100 ~ elig.peers.avg.100,
                cluster = "blkgrp",
                data = dat[elig.cs == 1 & post.rainwise == 0])

iv.corner <- feols(adopt.any * 100 ~ peers.100 + peers.corner.100 | blkgrp + year |
                     peer.adopt.any.100 + peer.any.corner.100 ~
                     elig.peers.avg.100 + corner.peers.avg.100,
                   cluster = "blkgrp",
                   data = dat[elig.cs == 1 & post.rainwise == 0])

iv.bus <- feols(adopt.any * 100 ~ peers.100 + peers.bus.100 | blkgrp + year |
                  peer.adopt.any.100 + peer.any.bus.100 ~
                  elig.peers.avg.100 + bus.peers.avg.100,
                cluster = "blkgrp",
                data = dat[elig.cs == 1 & post.rainwise == 0])

iv.all <- feols(adopt.any * 100 ~ peers.100 + peers.corner.100 + peers.bus.100 | blkgrp + year |
                  peer.adopt.any.100 + peer.any.corner.100 + peer.any.bus.100 ~
                  elig.peers.avg.100 + corner.peers.avg.100 + bus.peers.avg.100,
                cluster = "blkgrp",
                data = dat[elig.cs == 1 & post.rainwise == 0])

m.census <- feols(adopt.any * 100 ~ peers.100 | blkgrp + year |
                    peer.adopt.any.100 +
                    peer.adopt.any.100:hi.act +
                    peer.adopt.any.100:hi.inc +
                    peer.adopt.any.100:hi.own ~
                    elig.peers.avg.100 +
                    elig.peers.avg.100:hi.act +
                    elig.peers.avg.100:hi.inc +
                    elig.peers.avg.100:hi.own,
                  cluster = "blkgrp",
                  data = dat[elig.cs == 1 & post.rainwise == 0])

# Define variable labels --------------------------------------------------------
varnames <- c(
  "year" = "Year",
  "blkgrp" = "Block Group",
  "peer.adopt.any.100.ols" = "Peer Adoptions",
  "elig.peers.cs.100" = "\\# Eligible Peers (CS)",
  "elig.peers.100" = "\\# Eligible Peers (RG)",
  "peer.any.corner.100.ols" = "Peer Adoptions Corner",
  "peer.any.bus.100.ols" = "Peer Adoptions Bus",
  "peer.adopt.any.100" = "$\\widehat{\\text{Peer Adoptions}}$",
  "peer.any.corner.100" = "$\\widehat{\\text{Peer Adoptions Corner}}$",
  "peer.any.bus.100" = "$\\widehat{\\text{Peer Adoptions Bus}}$",
  "peer.adopt.any.100:hi.act" = "$\\widehat{\\text{Peer Adoptions}}$*High Active",
  "peer.adopt.any.100:hi.inc" = "$\\widehat{\\text{Peer Adoptions}}$*High Income",
  "peer.adopt.any.100:hi.own" = "$\\widehat{\\text{Peer Adoptions}}$*High Own"
)

# Preview Table A.3 -------------------------------------------------------------
etable(ivbase, iv.corner, iv.bus, iv.all, m.census,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       drop = c("peers.100"),
       fitstat = c("n", "ivfall", "sargan", "wh"))

# Save Table A.3 ---------------------------------------------------------------
style_tex <- style.tex(
  tablefoot.value = "",
  model.title = "",
  depvar.title = "",
  var.title = "\\midrule",
  stats.title = "\\midrule"
)

etable(ivbase, iv.corner, iv.bus, iv.all, m.census,
       tex = TRUE, digits = 3, depvar = FALSE, se.below = TRUE,
       dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       fixef_sizes = TRUE,
       drop = c("peers.100", "peers.bus.100", "peers.corner.100"),
       fitstat = c("n", "ivfall"),
       style.tex = style_tex,
       file = "output/tables/table_a3.tex", replace = TRUE)
